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ABSTRACT 

By constructing a global model based on 3D local magnctohydrodynamical (MHD) simulations, we 
show that the disk wind driven by magnetorotational instability (MRI) plays a significant role in the 
dispersal of the gas component of proto-planetary disks. Because the mass loss time scale by the 
MRI-drivcn disk winds is proportional to the local Kcplerian rotation period, a gas disk dynamically 
evaporates from the inner region with possibly creating a gradually expanding inner hole, while a 
sizable amount of the gas remains in the outer region. The disk wind is highly time-dependent with 
quasi-periodicity of several times Kcplerian rotation period at each radius, which will be observed as 
time-variability of protostar-protoplanetary disk systems. These features persistently hold even if a 
dead zone exists because the disk winds are driven from the surface regions where ionizing cosmic 
rays and high energy photons can penetrate. Moreover, the predicted inside-out clearing significantly 
suppresses the infall of boulders to a central star and the Type I migration of proto-planets which are 
favorable for the formation and survival of planets. 

Subject headings: accretion, accretion disks — MHD — stars: winds, outflows — planetary systems: 
formation — planetary systems: protoplanetary disks — turbulence 



1. INTRODUCTION 

Planets are believed to form in protoplanetary disks 
consisting of gas and dust components around newly 
born stars. The evolution of the gas component is cru- 
cial in determining the final state of the planetary sys- 
tem, such as the number and locations of terrestrial 
(rocky) and Jovian (gas-giant) planets (Ida & Lin 2004). 
The amount of the gas that is captured by Jovian plan- 
ets is generally much smaller than the total gas of the 
disk. Therefore, the gas component should dissipate 
via other mechanisms with the observationally inferred 
timescale of 10 6 — 10 7 yr (e.g., Haisch, Lada, & Lada 
2001; Hernandez et al. 2008). The currently favored sce- 
nario is that the gas dissipates by the combination of 
photocvaporation by Ultraviolet (UV) flux from a cen- 
tral star and viscous accretion (e.g. Shu, Johnstone, & 
Hollcnback 1992; Matsuyama, Johnstone, & Hartmenn 
2003; Takeuchi, Clarke, & Lin 2005; Alexander, Clarke, 
& Pringle 2006; Chiang & Murray-Clay 2007). How- 
ever, the time-evolution of the luminosity and spectrum 
of the UV radiation is quite uncertain, and moreover 
some observed transitional disks with inner holes are in- 
consistent with the photoevaporation mechanism (Calvet 
et al. 2005; Espaillat et al. 2008; Hughes et al. 2009). 

By performing local 3D ideal MHD simulations in the 
shearing box approximation (Hawlcy, Gammie, & Bal- 
bus 1995; Matsumoto & Tajima 1995), Suzuki & In- 
utsuka (2009, SI09 hereafter) have shown that MHD 
turbulence in protoplanetary disks effectively drives disk 
winds. The local timescale of the dynamical evapora- 
tion by magnetorotational instability (MRI) driven disk 
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winds is defined as 



£/(/w z )v 



(1) 



where (pv z ) w is the sum of the mass fluxes from the upper 
and lower disk surfaces and E is the surface density. We 
can estimate t cv ~ several thousand years at 1 AU of the 
typical protosolar disk (e.g., minimum mass solar nebula, 
or MMSN in short hereafter; Hayashi 1981); here we use 
the typical values, S = 1700 g cm -2 and (pv z ) w ~ 10~ 8 
g cm~ 2 s _1 (see §2.2 for the scaling of the disk wind flux). 
While the actual dispersal time is much longer as shown 
in this paper because the radial mass accretion is not 
taken into account here, this estimate shows that the 
MRI-driven disk wind should play a significant role in 
the dispersal of protoplanetary disks. In this paper, we 
investigate the evolution of gas density with disk winds 
and radial mass accretion in a global model. 

In SI09, we did not take into account the effects of 
a so-called dead zone which is inactive with respect to 
MRI because of the insufficient ionization for the cou- 
pling between gas and magnetic fields. The dead zone is 
believed to form around the midplane in the inner parts 
of protoplanetary disks because ionizing cosmic rays and 
X-rays cannot penetrate from the disk surfaces (Sano et 
al.2000). In this paper, we improve the models of SI09 
by investigating the effects of dead zones by performing 
resistive MHD simulations. 

The construction of the paper is as follows : We firstly 
describe the MHD simulations in local shearing boxes 
with and without dead zones in §2. In §3 we describe our 
global model. In §4, we show how the MRI disk winds 
disperse the gas component of protoplanetary disks from 
the inner part. Also, we show its effects on the planet 
formation. In §5, we discuss the properties of disk evo- 
lution mainly focusing on the escape of the disk winds 
from the gravity of a central star. 

2. LOCAL SIMULATIONS 
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Fig. 1. — Dynamical evaporation of a protoplanetary gas disk by local 3D ideal MHD simulation without mass supply by radial accretion. 
We impose weak vertical magnetic field, B z , with the plasma /3 value, /3 3!m jd = ^Pmid^/B^ = 10 6 , at the midplane. The lower right panel 
shows the total mass normalized by the initial mass in the local simulation box as a function of rotation time. The four color panels show 
the snapshots of the local protoplanetary disk simulation at t — (initial condition), 200 rotations, 2200 rotations, and 4130 rotations. 
The x, y, and z components respectively correspond to radial, azimuthal, and vertical components. The unit of each component is scale 
height, H = ^/2c s /n. The box size is (x,y,z) = (±0.5//, ±2H, ±4H ), which is resolved by (32,64,256) grid points. The colors indicate 
density normalized by the initial value at the midplane, the white solid lines illustrate magnetic field lines, and the arrows show velocity 
field. A small number of magnetic field lines (vertical lines) in the panel of t = reflect that the initially imposed magnetic field is weak 
(the number of field lines is scaled by magnetic field strength). We should emphasize that the actual dispersal is much slower because of 
the mass supply by accretion (see text). 



The main purpose of this section is to determine the 
turbulent viscosity and the mass flux of the disk winds 
which are used in the global models in §3. For that pur- 
pose, we perform 3D MHD simulations of a local proto- 
planetary (accretion) disk in the shearing box coordinate 
(Hawley, Gammie, & Balbus 1995). We use the simula- 
tion box with (x, y, z) = (±0.5H, ±2H, ±4ff), where the 
x, y, and z components respectively correspond to radial, 
azimuthal, and vertical components and scale height, H , 
is defined from sound speed, c s , and disk rotation fre- 
quency, fi, as H = \/2c s /il. 

In SI09 we have already shown results of ideal MHD 
simulation in the shearing box up to 400 rotation time. 
In this paper, we first extend this simulation until 5000 
rotation time when the significant mass is lost from the 
simulation box (§2.2). In §2.3 we take into account the 



effects of resistivity (i.e. dead zone). Later in this paper, 
we perform simulations in boxes with larger vertical ex- 
tents to study the effects of the box size on the escaping 
mass (§5.4.1). 

2.1. Launching of Disk Winds 

Before showing results of the local simulations, we dis- 
cuss basic properties of the disk wind obtained in SI09 
because the disk wind is a key that controls the evo- 
lution of protoplanetary disks in this paper. In SI09 
we interpreted that the disk winds are driven by the 
breakups of channel-mode flows (e.g. Sano et al. 2004) 
as a result of MRI (Balbus & Hawley 1991). Large- 
scale channel flows (e.g. Sano et al. 2004) develop most 
effectively at 1.5 - 2 times the scale height above the 
midplane. The breakups of these channel flows by mag- 
netic reconncctions drive disk winds in a time-dependent 
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manner with quasi-periodic cycles of 5-10 rotation pe- 
riod. The disk material itself is lift up recurrently, 
which will be observed as the time variation of effective 
disk surfaces. Time-variabilities are actually observed in 
protostar-protoplanctary disk systems (e.g. Wisniewscki 
et al. 2009; Muzerolle et al. 2009; Bary, Leisenring, & 
Skruskie 2009), which might be explained by quasi- 
periodic breakups of channel flows. The quasi-periodic 
feature of the disk winds is universally found not only 
in ideal MHD simulations but in simulations with dead 
zones as will be shown in §2.3. 

We should note that upward motions in vertically 
stratified accretion disks have been widely discussed with 
various interpretations. Magnetic buoyancy (Parker in- 
stability; Parker 1966) is one of the mechanism in up- 
lifting mass and magnetic field in the upper regions with 
\z\ > 1.5H (Miller & Stone 2000; Machida, Hayashi, 
& Matsumoto 2000; Nishikori, Machida, & Matsumoto 
2006). In fact, we observe -^-shape magnetic field struc- 
tures, which are characteristic of Parker instability, in 
our simulations as well (Figure 1). Recently, Latter, Fro- 
mang, & Gressel (2010) also discussed that channel flows 
are moving upward by magnetic buoyancy in a strati- 
fied disk. Magnetic buoyancy may play a role in upward 
flows and magnetic fields in cooperating with MPJ chan- 
nel flows in the upper regions (|z| > 1.5H). Detailed 
analysis of the role of magnetic buoyancy is driving the 
disk winds remains to be done. 

Flows and magnetic motions around the midplanc 
show complicated behaviors. On one hand, magnetic 
pattern seems to propagate away from the midplane 
(Davis, Stone, & Pessah 2010; Gressel 2010, ; our simu- 
lation shows the same trend.). Since this region is buoy- 
antly stable (Shi, Krolik, & Hirose 2010), this apparent 
propagating pattern may be due to a different mechanism 
from magnetic buoyancy. Brandenburg et al. (1995, see 
also Gressel 2010) try to interpret this upgoing pattern 
from dynamo waves (Parker 1955; Yoshimura 1975; 
Vishniac & Brandenburg 1997). On the other hand, 
the direction of the Poynting flux associated with mag- 
netic tension is toward the midplanc as we have shown 
in SI09. We explained that the breakups of large-scale 
channel flows at injection regions, z sa ±1.5if, drive mass 
motions to a midplane as well as surfaces. 

Although properties of vertical flows and magnetic mo- 
tions are not clearly understood in detail, upward flows 
in z > a few scale heights seem to be robust phenomena 
in stratified disks. The mass loading at the injection re- 
gions appear to be operated by the breakups of channel 
flows of MRI, and magnetic buoyancy also plays a role 
in lifting up the gas in the upper regions. 

2.2. Disks without Dead Zone 

We show results of the ideal MHD simulations in this 
subsection. We adopt the same resolution of the grid 
points, (x,y,z) = (32,64,256), as in SI09 for runs with 
different net vertical magnetic flux. In addition to the 
standard resolution runs, we perform the simulations 
with higher resolution, (x, y, z) = (64, 128, 512), in some 
cases. Extending from SI09, we carry out a long-time 
simulation in this paper until the significant mass is lost 
from the local box by the disk wind. Figure 1 is the result 
of the initial plasma j3 value, /3 2 . m id = 8-7T p m idCg / B 2 Z = 
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Fig. 2. — Dependences of averaged turbulent viscosity, a, (top) 
and mass flux of disk winds, C w , (bottom) on the initial plasma 
Az,mid(= 87rp m idCg/S^) values for the net vertical field, B z . The 
left-most grid corresponds to the initial toroidal field cases. The 
open symbols are the results of the standard resolution simulations 
with mesh points (x, y, z) = (32, 64, 256) and the filled symbols are 
the results of the higher resolution of (64, 128, 512). The solid lines 
represent <Sg and C w fl, which we use for the global model. The 

dashed lines denote the linear dependence on 



10 6 , at the midplane for net vertical magnetic field 4 , B z , 
where p m id is the density at the midplane. The figure 
shows that about 90% of the initial gas is dispersed in 
5000 rotations (5000 years at 1 AU) if there is no mass 
supply by accretion. This result shows that the evap- 
oration timescale using Equation (1) gives a reasonable 
estimate and the MRI-drivcn disk winds should signifi- 
cantly affect the evolution of protoplanetary disks. 

From the results of the local MHD simulations, we can 
determine the Shakura & Sunyaev (1973) a viscosity and 
the mass flux of the disk winds, a is calculated from 

4 Note that the integrated vertical magnetic flux, J dxdyB z is 
a strictly conserved quantity in the local shearing box simulations 
(Hawlcy, Gammie, & Balbus 1995), though the magnetic energy 
of z component f dxdyB^ /4-7T does not conserve. Then, , m id is 
a good indicator of magnetic field strength not only at the initial 
state but at later times when the magnetic field is amplified by 
MRI. 
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anisotropic stress of MHD turbulence, 

a = (v x 5v y -^L)/cl (2) 

where 5v y = v y + (3/2)f2x is the velocity shift from the 
background Kepler rotation (3/2)f2x). We use the aver- 
age a in the entire simulation box, 

_ _ J padxdydz 
J pdxdydz 

Note that the density weighted average, a, is directly 
related to mass accretion rate (e.g., Pringle 1981). As for 
the disk wind mass flux, we use the nondimensionalized 
mass flux, 

C w = (pv z ) w /(p mid c s ). (4) 

Figure 2 shows a and C w as functions of /3 z ,mid- Note 
that the net B z becomes strong from left to right. The 
left-most grid is for the cases of no net B z field. In these 
cases we initially give purely toroidal field, B y , with the 
j3 values of 10 6 in — 3H < z < 3H, whereas the results 
do not depend on the initial strength because magnetic 
flux of y (as well as x) component does not conserve. 

The figure exhibits that both a and C w have the floor 
values for sufficiently weak net vertical magnetic field, 
As. mid 10 6 . In these cases the saturated magnetic field 
at the midplane gives (3(= 8irp/B 2 ) ~ 100, indicating « 
1% of the gas energy is transferred to the magnetic fields 
(SI09). The saturation level roughly corresponds to the 
floor value of a ~ 0.01 (a is roughly 1/13). Such weak net 
vertical field gives little effects, because it is much smaller 
than the turbulent component of the vertical fields. The 
turbulent component of B z gives {B 2 )/8ir ~ 10~ 4 — 10~ 3 p 
even in the zero net vertical field case, where () denotes 
time-average. This value seems to determined as the level 
that is one or two orders of magnitude smaller than the 
dominant toroidal component ((B 2 ) /8n ~ p/f3 ~ O.Olp). 
If the net vertical field is much smaller than the turbulent 
vertical component, (B z ) 2 /8n < 10~ 6 p the values of a 
and C w are not affected by the net vertical field. 

On the other hand, for larger magnetic field cases, 
As, mid ^ 10 5 , the net vertical magnetic field plays a 
role, because it is not negligible compared to the turbu- 
lent component, a and C w increase almost linearly with 
magnetic energy of net vertical field (cx l//3 2 , m id)- The 
behaviors of a and C w are similar, which indicates that 
the mass flux of the disk winds is positively correlated 
with the mass accretion rate. This is reasonable because 
the energy source of the disk winds is the gravitational 
energy liberated by accretion, which is discussed in §5. 

In the higher resolution runs, the saturation levels of 
the magnetic fields are slightly lower than the corre- 
sponding cases with lower resolutions. Then, both a and 
C w become slightly smaller in the higher resolution runs. 
The dependence of the saturation level of magnetic field 
and a on resolutions is still under debate and widely dis- 
cussed in various authors (e.g. Balbus & Hawley 1998; 
Pessah, Chan, & Psaltis 2007) 

2.3. Disks with Dead Zone 

The temperature of protoplanetary disks are too low 
to ionize the gas by thermal collisions. Various ionization 



sources, such as X-rays, cosmic rays, and radioactive nu- 
clei, have been widely discussed (Umebayashi & Nakano 
1980; Hayashi 1981; Glassgold, Najita, & Igea 1997). 
The ionization degree at midplanes is generally smaller 
than that in upper regions because the recombination 
rate is higher there as a result of higher density. Un- 
der certain circumstances, dead zones, in which the gas 
is decoupled with the magnetic fields due to the insuffi- 
cient ionization, are supposed to form near midplanes. 
However, Inutsuka & Sano (2005) also introduced a 
self-sustained mechanism by current-carrying electrons 
in turbulent disks, which increases the resultant ioniza- 
tion degree. Thus, calculating ionization degree is not 
straightforward. In this paper, we study extreme cases 
which give distinct dead zones. As will be described be- 
low, we take into account the ionization by cosmic rays 
and X-rays. Since they come from disk surfaces, dead 
zones tend to form near the midplane if the surface den- 
sity is sufficiently high. We adopt a simple treatment in 
determining ionization degree and study effects of dead 
zones on disk winds. 

2.3.1. Set-up 

To treat dead zones, we take into account resistivity 
in the induction equation that describes the evolution of 
magnetic fields : 

f) /■? 

— = V x (v x B -jjV x B), (5) 

where rj is resistivity which is determined by ionization 
degree, x c as n = 234VT/xe (Blaes & Balbus 1994), 
because electrons control the coupling between gas and 
magnetic field in high density medium. We assume the 
temperature structure of the MMSN (Hayashi 1981) : 

t = 293K (tau) < ^ 

which can be transformed into the sound speed, c s = 




-35 -30 -25 

log(f/n) 



Fig. 3. — Ionization degree, x e , on ionization rate per number 
density, f/n. The data are taken from Inutsuka & Sano (2005); 
Sano et al. (2000). 
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The ionization degree is determined by the balance 
between ionization and recombination. In this paper, 
we use the result of previous calculation by Sano ct al. 
(2000) and Inutsuka & Sano (2005), which give x e as a 
function of ionization rate, £, and local number density, 
n, by calculating recombination on dust grains and ra- 
diative and dissociative recombination in gas phase (see 
also, e.g. Umebayashi & Nakano 1980; Ilgner & Nelson 
2006). When one fixes the abundance of gas and the 
properties of dust grains, the ionization degree is a func- 
tion of the only one parameter, £/n (Inutsuka & Sano 
2005). This is because the recombination is essentially 
two-body reaction and recombination rate per volume is 
proportional to n 2 while the ionization rate per volume 
is £n. Figure 3 presents x e for the solar abundance gas 
with gas-to-dust ratio of 100 and dust grain size of 0.1/xm 
(Inutsuka & Sano 2005; Sano et al. 2000). We use this 
data for our local resistive MHD simulations. We adopt 
the mean molecular weight, /i = 2.3, for the conversion 
between n and p. This /i value reflects the condition that 
the major component is H2 molecules. 

We take into account the ionization by Galactic cosmic 
rays and X-rays from a central star. The total ionization 
rate, £, is the sum of these two sources, 

£ = £cr + £x, (7) 

where subscripts CR and X represent cosmic rays 
and stellar X-rays, respectively. We adopt £cr = 
10~ 17 exp(— l/l cr ) s _1 , where I is the column densities in- 
tegrated from the upper or lower surfaces and l cl - = 100 
g cm~ 2 is the path length of cosmic rays (Hayashi 1981; 
Turner, Sano, & Dziourkevitch 2007). 




Fig. 4. — The path, ds, of emitted X-ray from the source located 
above a central star to the disk (pink line). The blue line is the 
upper disk surface, (see text) 

Observations of T-Tauri stars show high X-ray ac- 
tivities. The typical luminosity in X-ray wavelength is 
10 29-31 erg s _1 , which is 3-5 orders of magnitude higher 
than the level of the present Sun (e.g. Telleschi et al. 
2007). These X-rays are supposed to be an efficient 
ionization source. Following Glassgold, Najita, & Igea 
(1997), we model the ionization rate, £x, by the X-rays. 
We assume the X-ray sources located at 10 Rq (~ 3 — 5 
stellar radii) above and below a star, where Rq denotes 
the solar radius. The X-ray follows the path shown in 
Figure 4. The ionization rate is obtained by Glassgold, 
Najita, & Igea (1997); Fromang, Terquem, & Balbus 
(2002) as follows: 




Fig. 5. — Snap-shot of the local disk with dead zone at 250 rota- 
tions in the case with medium X-ray activity at 1 AU. The colors 
indicate the logarithmic scale of density, the solid lines represent 
magnetic fields, and the arrows are velocity field. 
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where Lx is the X-ray luminosity, a ~ 4 x 
10~ 24 (jjfjy) is the absorption cross section depend- 
ing on the X-ray energy, Ex 7 t = J dsna is opti- 
cal depth and J(t) = 0.6867- - 606 exp(-1.778r a262 ). 
Here, the optical depth, r, is integrated along the path, 
ds, in Figure 4 for the variation of scale height, H w 

0.05AU( T i u ) 5/4 . 

We perform simulations at three different locations, 
r = 1,5,25 AU. We consider the three levels of X- 
ray activities : (X x (erg s" 1 ), E x QtcV)) = (10 31 ,5), 
(10 30 ,3), and (10 29 ,1), which we call the strong, 
medium, and weak X-ray cases, respectively. We 
adopt the original MMSN for the surface density, £ = 

1700 g cm~ 2 (jjjj) ~ 3/2 - Iu the resistive MHD simula- 
tions we adopt the only standard resolution ((x,y,z) = 
(±0.5H,±2H,±4H) is resolved by (32,64,256) grid 
points). We impose weak net vertical magnetic field, 
/3 z ,mid = 10 6 , at the midplane. The constant /3 2jm id dis- 
tribution indicates that we assume the dependence of 

the net vertical field of (B z ) m 0.01 G {jxu)~ 13/8 for the 
MMSN {p mid oc r- 11 / 4 and c s cx r" 1 / 4 ). 



6t/J(r) = 1-2 x lQ-^cm-V 1 ) 



2.3.2. Result 




Fig. 6. — Comparison of time-height diagram of pv z between the resistive MHD case of the medium X-ray (Lx = 10 30 erg s 1 and 
Ex = 3 keV) at 1 AU (left panel) and the ideal MHD case (right panel). pv z is averaged over the x — y planes. The horizontal axis is in 
unit of rotation period, and the vertical axis is normalized by scale height, H = \/2c s /0. 



First, we describe results of the case with the medium 
X-ray activity at 1 AU in detail before discussing results 
of the different locations and X-ray activities. Figure 5 
shows a snapshot of the local disk structure at 250 rota- 
tions. (Note that the density is in logarithmic scale here; 
c./., a linear scale was used in Figure 1.) The hgure shows 
that the magnetic field lines become almost straight in 
the dead zone, —2H < z < 2H, because MRI does not 
operate due to the insufficient ionization here. This is 
a clear contrast to disks without dead zones (Figure 1; 
sec also SI09). In the surface regions, however, the mag- 
netic fields are more turbulent because the ionizing cos- 
mic rays and X-rays can penetrate to these regions and 
MRI is active at r « ±(2 — 3)H. One can also see that 
the disk winds stream out of both surfaces because the 
disk winds are driven from the surface regions with suffi- 
cient ionization rather than a deeper midplane location. 
The breakups of the channel flows triggered by MRI at 
z ps 2H and Parker instability in z > 3H drive these disk 
winds, as explained in §2.1. For example, one can see a 
typical D-shape channel flow structure at z ~ —2H, and 
a ^-shapc structure at z ~ (3 — 4) if, which is typical 
for Parker instability. Although the mass flux of the disk 
winds become moderately smaller than the ideal MHD 
run, the dead zone give little effects on the disk winds 
(see below). 

Figure 6 compares the time-height diagrams of the 
mass flux, pv z of the same case (left panel) to the result 
of the ideal MHD case (right panel) . pv z is averaged over 
the x — y planes at each time step. The left panel shows 
that disk winds flow out of the upper and lower sur- 
faces even though the dead zone forms around the mid- 
plane. This is because the disk winds are excited from 
the surface regions, which we called injection regions in 
SI09, with the sufficient ionization, rather than deeper 
locations near the midplane. These injection regions 



are a consequence of the breakups of large-scale chan- 
nel flows. The disk winds flow out intermittently with 
quasi-periodicity of 5-10 rotations owing to the quasi- 
periodic breakups of the channel flows. The recurrent 
nature of the disk winds universally holds in both cases 
with and without dead zones, which might explain ob- 
served time- variabilities of protostar-protoplanetary disk 
systems (Wisnicwscki et al. 2009; Muzerolle et al. 2009; 
Bary, Leisenring, & Skruskie 2009). The heights of the 
injection regions become slightly higher than those of 
the ideal MHD case because the ionization degree at the 
deeper locations is not sufficient for MRI. Accordingly 
the mass flux of disk winds from the dead zone case is a 
little smaller than that of the no dead zone case. 

Figure 7 compares the vertical structure of these cases 
averaging over 200 rotations after the quasi-steady states 
are achieved. a z in the lower right panel is the average 
of a (Equation 2) on each x — y plane. The right panels 
(magnetic energy and a values) show that the dead zone 
extends from z = —2H to 2H in the resistive MHD case. 
a z sharply declines at z = ±2H from the surface regions 
toward the midplane. It is expected that the mass accre- 
tion mainly takes place in the surface regions (Gammic 
1996). The averaged a value, a(= J dzpa z /^,), which di- 
rectly determines global mass accretion rate, of the dead 
zone case is ~ 3 x 10 -4 , while it is ps 0.011 in the no dead 
zone case (Figure 2). On the other hand, in the surface 
regions the magnetic energy and a z of the resistive MHD 
case are similar to those of the ideal MHD case because 
the sufficient ionization is achieved there by the X-rays 
and cosmic rays from the surfaces. The disk winds are 
effectively accelerated from the injection regions, and the 
velocity structures are very similar in both cases (the up- 
per left panel of Figure 7). The density of the winds are 
smaller because the heights of the injection regions are 
slightly higher in the dead zone case. Accordingly, the 
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Fig. 7. — Comparison of the resistive MHD (dead zone) case 
of the medium X-rays at 1AU (solid lines) with the ideal MHD 
(no dead zone) case (dashed lines). Both cases assume the net 
vertical magnetic field with /? 2) mid — 10 6 . The horizontal axis 
is vertical height in unit of H(= \/2c s /lQ). The upper left panel 
shows the vertical velocities normalized by the sound speed; the 
lower left panel shows the densities; the upper right panel shows 
the magnetic energy. The lower right panel shows a z and pa z . 



mass flux of disk winds of the resistive MHD case is the 
half of the mass flux of the ideal MHD case. 

Figure 8 summarizes the averaged turbulent viscosities, 
a, and the mass fluxes of the disk winds for the different 
locations and X-ray activities. The figure illustrates that 
both a and C w are smaller in the inner regions. Since 
surface density is larger for smaller r, ionizing cosmic 
rays and X-rays cannot reach the midplane in the inner 
disk. A large fraction is occupied by the dead zone in the 
inner parts of disks. The a values are reduced by more 
than an order of magnitude at r = 1 AU. On the other 
hand, C w is not so reduced, being half at most at r = 1 
AU, because the disk winds are excited from the surface 
regions and not so affected by the dead zones. 



3. GLOBAL MODELING 

We solve the time-evolution of surface density with 
mass accretion and disk wind mass loss (e.g. Alexander, 
Clarke, & Pringle 2006): 



~dt 



ld_ 

r dr 



2 d , 2 



+ (pv z ) w = 0. (9) 
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Fig. 8. — Turbulent viscosity, a, (top) and disk wind mass flux, 
C(= (pf z )w/(Pmid c s)), (bottom) of the resistive MHD simulations 
at different locations. The squares are the results of the strong X- 
ray cases (Lx = 10 31 erg s — 1 , E x = 5 keV), the triangles are the 
results of the medium X-ray cases (L x = 10 30 erg s _1 and E x = 3 
keV), and the circles are the results of the weak X-ray cases (Lx = 
10 29 erg s — 1 and Ex = 1 keV). The solid lines are the fitting 
formula (Equations 14 and 16) used for the global calculations. 
The dotted lines arc the results of ideal MHD simulations with net 



vertical field, fi z „ 



10" 



The second term denotes the radial mass flow by turbu- 
lent viscosity. Here and from now, we simply write a 
for the turbulent viscosity in the global models, which is 
adopted from a of the shearing box simulations of the 
previous section. The third term is the mass loss by disk 
winds, where we here assume the specific angular mo- 
mentum carried in the disk winds is the same as that in 
the disk material. Both a and (pv z ) w are adopted from 
the local 3D MHD simulations in §2.2 and 2.3. 

We assume the temperature structure of the MMSN 
(Equation 6). The initial surface density profile is also 
adopted from the MMSN, 



/ r \ -3/2 
S = ^ sSo VTAU ) e M~r/r. 



(10) 



where S = 2400 g cm" 2 at 1 AU (/ g = 0.7 for the 
original Hayashi MMSN) and we use a cut-off radius, 
r C ut = 50 AU. Although in the local resistive MHD sim- 
ulations for dead zones (§2.3) we adopted the original 
MMSN value, f g ^o = 1700 g cm" 2 , a specific choice of 
f g in the local simulations does not change the results of 
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the global disk because we explicitly take into account 
the dependences of a and C w on E (§3.2). 

We integrate Equation (9) by using the nondimension- 
alized variables in unit of ro = fio = ZgEo = 1, so the 
scaling factor, / g , does not appear explicitly in the calcu- 
lations. The calculation region 5 is from r ln = 0.01 AU to 
fout = 10000 AU which is resolved by 2000 mesh points 
with grid spacing in proportion to ^Jr. 

3.1. Disks without Dead Zone 

We apply the a and C w obtained in the local ideal 
MHD simulations (§2.2 and Figure 2) to the global 
model. Figure 2 can be regarded as a time sequence along 
with disk evolution because the surface density of pro- 
toplanctary disks decreases while vertical magnetic flux 
is supposed to be rather kept constant. If initial vertical 
magnetic fields arc zero or very weak (/3 Zi mid 3> 10 7 ), 
a and C w are expected to stay almost constant even 
after surface density decreases considerably. For such 
situations, we use constant a = an = 8 x 10~ 3 and 
C w = C Wi fi = 2 x 10~ 5 (solid lines in Figure 2), which 
are the floor values obtained by the local simulations. 
For the disk wind mass flux, we choose a conservative 
value because the actual mass flux might be moderately 
smaller by returning mass from higher altitudes (§5.1 and 
5.4). 

When C w stays constant, the disk wind flux has the 
following scaling : 

{pv z ) w = C w p mid c s cx Er~ 3/2 , (11) 

where for the last proportionality we use p m idC s oc Ef2 
and assume a Kcplcrian rotating disk . Equation (11) 
shows that the wind mass flux is larger for smaller r, and 
the dispersal of protoplanetary disks by the disk winds 
starts from the inner part. 

If initial vertical fields are not so weak, a and C w even- 
tually increase when /3 z , m id 10 5 (Figure 2). In order 
to take into account this effect we adopt the following 
prescription : 

a = a R x max(l, %^), (12) 
E(r) 

and 

C w = C w . fl x max(l, %^), (13) 
E(r) 

where E up is the surface density at which a and C w start 
to increase in Figure 2 {fi z m \i ~ 10 5 ). E up is deter- 
mined by the initial vertical magnetic flux. We model 
^up( r ) — ^upEi n it(f); oi and C w starts to increase when 
the surface density decreases to <5 up of the initial value. 
For simplicity, we assume a constant <5 up = 0.01 in this 
paper. 

We calculate the three cases for the disks without dead 
zones summarized in Table 1 (Models I - III). In Models 

5 Our results are not affected by the location of r- ln . The adopted 
value, r; n « 2Rq, roughly coincides a typical radius of T-tauri 
stars, where Rq is the solar radius. In realistic situations, how- 
ever, a disk may truncate at several stellar radii where the Kep- 
lerian rotation frequency equals to the rotation frequency of the 
corotating stellar magnetic field (Ghosh & Lamb 1979). The mat- 
ter can directly accrete to a central star through connecting flux 
tubes (Kenyon, Yi, & Hartmann 1996). 



Model 


Disk Wind 


net B z 


Dead Zone 


I 


No 


Weak/No 


No 


II 


Yes 


Weak/No 


No 


III 


Yes 


Strong 


No 


IV 


No 


Weak/No 


Yes 


V 


Yes 


Weak/No 


Yes 



TABLE 1 
Global models. 



II and III we take into account the disk winds; Model 
II adopts the constant a & C w to model weak vertical 
magnetic fields and Model III prescribes Equations (12) 
and (13) to incorporate relatively strong vertical fields. 

3.2. Disks with Dead Zones 

We apply the results of the local resistive MHD sim- 
ulations (§2.3 and Figure 8) to the global model. An 
essential point is that a and C w can be determined by r 
and E. Then, we use the following parameterization : 

a = a s J^- 1 (14) 

where E act is the sum of the column density of active 
regions near upper and lower disk surfaces which is mod- 
eled as 

E act =minfs CR + E x (^)~ 2 ,EY (15) 

where Ecr is the column density with sufficient ioniza- 
tion by cosmic rays, and Ex is the column density with 
sufficient ionization by X-rays normalized at 1 AU. Ex 
has dependence on r~ 2 to take into account the geo- 
metrical dilution of the X-ray flux from a central star, 
while Ecr is constant because cosmic rays are diffusely 
distributed. To reproduce the results of the local sim- 
ulations, we adopt Ecr = 12 g cm -2 and Ex = 25 g 
cm~ 2 . The mass flux of disk winds is also reduced in 
accordance with a, but we set a lower limit to match the 
local simulation results: 

(X 

Cw Cw,min (Cw.fl Cw,min) (1^) 
where We USe C w m in 

= 0.45 x C w fl. 
The solid lines in Figure 8 represent Equations (14) and 
(16). In the figure we use a n = 0.011 and C w £ — 7.7 x 
10 -5 for the floor values of the ideal MHD simulations 
with net vertical field of /3 z . m id = 10 6 (Figure 2). For 
the global disk calculation we use the original values, 
a a = 8x lO- 3 and C w , fl = 2 x 10" 5 . 

4. RESULTS OF GLOBAL MODELS 

4.1. Evolution of Gas Disks 

The top panel of Figure 9 shows the evolution of the 
surface densities of the no dead zone cases (Models I - 
III). The result of the no wind case (Model I; black thin 
lines) follows a self-similar solution of E oc 1/r in the 
inner region with an exponential cutoff in the outer re- 
gion (Lyndcn-Bcll & Pringlc 1974). On the other hand, 
the disk wind cases (Models II and III; blue thin and 
red thick lines) shows faster decreases of S in the inner 
regions owing to the contribution from the disk winds in 
addition to the accretion. The figure clearly shows that 
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r(AU) 



Fig. 9. — The results of the no dead zone cases (Models I - III). 
From top to bottom, time evolution of disk surface density, E, 
timescale, Tdrifti °f inward drift of a m-size boulder, and timescalc, 
T mi g ,l; °f tyP e I migration of an Earth-mass planet are shown. 
The black thin lines are the results without disk winds (Model I); 
the blue thin lines are the results with disk winds for weak/no 
vertical magnetic fields (Model II); the red thick lines the results 
with disk winds for relatively strong net vertical fields (Model III). 
The dotted lines are the initial values. The dashed, solid, and dot- 
dashed lines are the results at 10 s , 10 6 , and 10 7 yrs. Note that E 
and T m ig,l can be scaled by / g . For example, the case with / g = 2 
gives twice larger S and smaller T m ; g j than the case with / g = 1. 
The results of 7d r ;f t are for / g = 1, because this scaling cannot be 
applied to Thrift- The bending points of Thrift correspond to the 
change of the regimes of drag force. The inside region corresponds 
to the Stoked regime where the dust size is larger than the mean 
free path of a gas particle, and the outside region corresponds to 
the Epstein regime where the dust size is smaller that the mean 
free path of a gas particle. 

the dispersal of the gas disk takes place in an inside-out 
manner, as discussed above (Equation 11). 

The case with relatively strong vertical net magnetic 
flux (red lines) shows an expanding inner hole because 
the surface density decreases to reach 5 up (= 0.01) of the 
initial values faster in the inner regions and a and C w 
start to increase from inner to outer locations. Although 
the quantitative properties of evolving inner holes de- 
pends on the adopted parameters (see §2.2), the observed 
properties of transitional disks with inner holes (Calvet 
et al. 2005; Espaillat et al. 2008; Hughes et al. 2009) may 
be explained by the mechanism presented here. 

The top panel of Figure 10 shows the evolution of the 
surface densities of the dead zone cases. The results of 
Models IV & V exhibit density enhancements, which cor- 
respond to the dead zones. Because the surface densities 




O,- -, 

0.01 0.1 1 10 100 1000 

r(AU) 

Fig. 10. — Comparison of the dead zone cases with the no-dead 
zone case. The black, red, and blue lines are the results of Models 
IV (dead zone / no disk wind), V (dead zone / disk wind), and II 
(no dead zone / disk wind), respectively. The top panel compares 
the time evolution of surface density, E. The dotted, dashed, solid, 
and dot-dashed lines are the results at t = 0, 10 s , 10 6 , and 10 7 
yrs. The middle panel compares the inward drift timescales, Thrift, 
of a m-size boulder at t = 10 6 yr. The bottom panel compares the 
timescale, r m i g: i, of type I migration of an Earth- mass planet at 
t = 10 6 yr. 

in these regions are high, the X-rays and cosmic rays can- 
not penetrate to the midplanes. The dead zones form 
around the midplanes, and a and C w become smaller. 
Smaller a leads to slower mass accretion 6 , and then, the 
mass accumulates around the dead zones, which is ob- 
served as the density enhancements. Without disk winds, 
the dead zone exists until 10 7 yr (Model IV; black lines). 

On the other hand, when taking into account the disk 
winds, the dead zone almost disappears at 10 6 yr (Model 
V; red lines) because the surface density decreases faster 
owing to the disk winds, which further leads to the effec- 
tive penetration of the ionizing X-rays and cosmic rays 
to the midplane. After the dead zone disappears, the 
disk evolution follows the no dead zone case. The sur- 
face density structure in later times (e.g. at 10 7 yr in the 
figure) is very similar to that of the no dead zone case 
(Model II; blue lines). We can conclude that the MRI- 
driven disk winds play an essential role in the dispersal 
of protoplanetary gas disks even though the dead zone 
forms. 

6 The disk expands around the outer edge of the dead zone, while 
the mass accretes in the rest of the dead zone region. 
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r(AU) 



Fig. 11. — Upper: Surface density of dust (solid) and gas (dashed) 
components at t = 10 6 yr. Lower. Dust-to-gas ratios at t = 10 5 yr 
(dashed), 10 6 yr (solid), and 10 7 yr (dot-dashed). The dotted line 
is the initial condition (dust-to-gas ratio = 0.01). 

4.2. Dynamics of Boulders and Planetesimal Formation 

In addition to the dynamical evaporation of proto- 
planetary gas disks, MRI-driven disk winds affect the 
planet formation at various stages. At an early stage the 
rapid infall of boulders to a central star (Weidenschilling 
1977), which hinders the growth of solids to larger bodies 
by aggregation, is a severe problem. The solid compo- 
nent rotates with the Keplerian velocity as a result of the 
force balance between the gravity due to a central star 
and the centrifugal force. On the other hand, the gas 
rotates with sub-Keplerian velocity by the contribution 
from outward pressure gradient force. Then, the rotation 
of the solid component is slightly slowed down because 
of the head-wind from the gas, and drifts inward. Under 
the typical MMSN condition, the infalling timescale of 
~ meter-sized boulders is ~ 100 — 1000 years at 1 AU, 
which is too rapid to form planetesimals (~ kilometer 
size) in a turbulent gas disk. 

However, our calculations show that the surface den- 
sity is increasing with r in the inner region (the top pan- 
els of Figures 9 and 10) . The inward drift rate becomes 
smaller than previously discussed. The middle panels of 
Figures 9 and 10 compares the inward drift timcscalcs 
(Weidenschilling 1977), 

/ i dp t s n y 1 
Tdnft = r ^*TTwj ' {) 

where t s is stopping time of solid material in a gas disk. 
Here we consider a one meter-size spherical boulder for 
t s , and the pressure-gradient force (dp/dr < for sub- 



Keplerian rotating gas disks) is estimated at the mid- 
planes of disks. As expected, Thrift becomes longer when 
the disk wind is taken into account (Figure 9), which is 
more favorable for the formation of planetesimals. The 
tendency is the same for the dead zone cases as well (Fig- 
ure 10), while Tdrift's show complicated behaviors at the 
edges of the dead zones reflecting the sharp density gra- 
dients. 

The dispersal of the gas component directly leads to 
the increase of a dust-to-gas ratio, which is also impor- 
tant in the context of gravitational instability of dust 
particles. Sekiya (1998) investigated the turbulence due 
to shear motions between gas and dust. He found that 
when a significant fraction of the gas component is dis- 
persed with dusts left, the shear-induced turbulence is re- 
duced and dusts become gravitationally unstable, which 
possibly leads to the formation of planetesimals (see also 
Youdin & Shu 2002; Michikoshi & Inutsuka 2006, for 
related works). Johansen et al. (2007) also proposed 
that streaming instability triggers the direct formation 
of large planetesimals or dwarf planets when the dust- 
to-gas ratio increases to an order of unity. The disk 
winds disperse the gas component selectively and raise 
the dust-to-gas ratio in an inside-out manner, which ac- 
tivates streaming instability from inner regions. 

To illustrate the increase of the dust-to-gas ratio, we 
calculate the time-evolution of dust surface density by 
using the same global model. Here, we assume that dust 
particles follow the only viscous accretion with gas near 
the midplane and do not escape with disk winds; we 
solve Equation (9) by setting (/w z ) w = for the dust 
component. This assumption is reasonable if dust par- 
ticles are well-coupled with gas near the midplane and 
decoupled in the upper regions. Nondimensional stop- 
ping time, flt s , is a good indicator which measures the 
coupling between dust and gas; if ttt s < 1, dusts are well- 
coupled with gas, and vice versa. When one takes dust 
particles with size of millimeter at 1 AU of the MMSN 
with f g = 1 (Equation 10) as an example, flt s «2x 10 -4 
at the midplane. Since stopping time is inversely pro- 
portional to gas density, Qt s exceeds unity at z ~ ±3_ff 
where pj 'pmid ~ 2 x 10~ 4 (Figure 7). The assumption is 
reasonable for moderately small dusts (sub-millimeter - 
meter at 1 AU of the MMSN). Further smaller dusts will 
be lift up with disk winds, while larger solid particles are 
decoupled with gas even at the midplane. 

We initially impose a uniform dust-to-gas ratio = 0.01 
and follow the evolution of both gas and dust with using 
the parameters of Model II of Table 1. The upper panel 
of Figure 11 presents the surface densities of dust (solid) 
and gas (dashed) at t = 10 6 yr, which shows that the dust 
surface density is larger than the gas surface density in 
the inner region, r < 0.1 AU. The lower panel shows that 
the dust-to-gas ratio gradually increases from the inner 
region. 

4.3. Type I Migration 

After the formation of planets, the gas component in a 
disk also plays a role in the evolution of the planetary sys- 
tem. A lower-mass planet, which cannot create a gap in 
a gas disk, resonantly interacts with the gas component 
of a disk through gravitational torque (Ward 1997). As 
a result planets generally migrate inward with timescale 
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Fig. 12. — Orbital evolution of planets by type I migration. The 
horizontal axis shows radial distance and the vertical axis shows 
time elapsed from the starting point of the disk calculation. The 
solid lines are the results with disk winds (Model II) and the dashed 
lines are the results without disk wind (Model I). The circles in- 
dicate the initial locations of newly formed planets : At t = 10 5 
yr and 10 6 yr we put planets with 0.3Mj at 0.3 AU, planets with 
lMfg at 1 AU, and planets with 5Mfg at 5 AU. 

(Tanaka, Takeuchi, & Ward 2002), 



t (r) « 5 x 10 V 



4.35 \ /E(r)\ _1 / M 

2.7+ 1.1s J \ Me 

(18) 

where s is the local gradient of S oc r s , M is planet 
mass, Mq is the Earth mass, and we assume the MMSN 
around a central star with the solar mass. This equation 
indicates that the migration is faster in a more massive 
(larger S) gas disk because of the larger torque on a 
planet. In typical situations, r m i g) i is shorter than the 
lifetimes of protoplanetary gas disks; newly formed ter- 
restrial planets and cores of gas-giant planets quickly fall 
into a central star. 

However, we have shown that the gas density in the 
inner regions quickly decreases by the disk winds, and 
the gradient of surface density becomes positive in the 
inner region. Consequently, T m i gi i becomes considerably 
longer than previously considered. The bottom panels of 
Figures 9 and 10 show r m i g j for an Earth-mass planet 
around a solar-mass star. The figures illustrate that the 
disk winds greatly suppress the migrations in the inner 
region both in no dead zone and dead zone cases. Simi- 
larly to Tdrift, T m i gi i's of the dead zone cases show com- 
plicated behaviors at the edges of the dead zones (Mat- 
sumura, Pudritz,& Thommes 2007; Kato et al. 2009). 

To illustrate the suppression of type I migration more 
clearly, we calculate the migrations of planets with the 
evolution of protoplanetary disks by using migration 
speed, r/r m i g j (Figure 12). In the no disk wind case 
(Model I) all the planets infall to a central star. On the 



other hand, in the disk wind case (Model II) the migra- 
tions are slow and all the planets survive. When the disk 
winds are considered, type I migration becomes unim- 



portant especially at later times, t > 10 5 — 10 6 
the typical MMSN condition. 



yr 



under 



5. DISCUSSIONS 

5.1. Energetics 

In this paper we have applied results of the local simu- 
lations to the global models. We use the mass flux of the 
disk winds at the upper and lower boundaries, z = ±47?, 
of the simulation box. Since the wind velocities at the 
boundaries are still smaller than the escape speed from 
a central star, we should carefully examine whether the 
disk winds really escape from the disks. In this subsec- 
tion, we examine the energetics of the disk winds to see 
whether the accretion energy can potentially drive the 
disk winds that can escape from the gravity of a central 
star. 

The starting point here is an equation for the conser- 
vation of total energy : 



d_ 

dt 



i 



+v- 



i 



p 



GAL 



7- lp 



pv | -v 



£2 

7 P GNh 



7- lp 



B 2 , 
+v — - —(B ■ v] 



-qio 



(19) 



where q\ oss is energy loss which is modeled below and here 
we consider the only r derivative. We integrate Equation 
(19) with J dz by neglecting the small terms except the 

anisotropic stress, Sv,pv r — B^Br/Airp = ac 2 s , whereas we 
separate into Keplerian rotation plus small pertur- 
bation, V0 = rfl + Svrf,. We also assume that accretion 

velocity and sound speed (c s = \pjpj~p) are smaller than 
rotation velocity, v r ,c s rQ. Then, the total energy of 
a ring at r changes according to 



d_ 

dt 



+- 



ld_ 

dr 



d 

rn—(r 2 Y,ac 2 s ) + r 2 Y>ttac 2 s 
dr 



= — <^lo 

(20) 

where Qioss = J q\ oss dz and we have used = rQ 2 

and r£iv = ■^■§^(r 2 Y,ac 2 ). The spatial derivative term 
on the left-hand side represents the energy liberated by 
mass accretion and viscous heating. 

We investigate whether this liberated gravitational en- 
ergy is sufficient to drive disk winds. For simplicity, we 
only consider the kinetic energy of disk winds for the loss 
term and neglect additional effects such as heating by 
UV/X-ray radiation (energy input), acceleration by stel- 
lar winds (momentum input), and radiative cooling (en- 
ergy loss) . Since we consider the disk winds from a Keple- 
rian rotating disk, we can write <5i OS s = \pv z {v 2 + r 2 Q 2 ). 
Disk winds can escape to infinity if v z exceeds the escape 
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speed, t>osc — V%rfl; namely if the condition, 
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r r 2 n 2 ' 


1 8 ' 


dt 


S 2 


r dr 



rft — (r 2 £ac 2 ) + r^Qac 2 
dr 



-\{pv z )^r 2 tf > (21) 

is satisfied, disk winds can be driven solely by the liber- 
ated gravitational energy of accretion, because the left- 
hand side of Equation (21) is the energy flux of disk winds 
at infinity. The time derivative term in Equation (21) can 
be eliminated by using Equation (9). Then, using the re- 
lation of Keplerian rotation, D, oc r -3 / 2 , Equation (21) is 
rewritten as 

^fffiac 2 > 2r 2 n 2 {pv z ) w , (22) 
or more specificly, 

8 ci( ro n y _ 

a 2 cU r ° = r *" (23) 

where we are using Equation (4) and the relation of the 
MMSN (c s oc r" 1 / 4 ) with Keplerian rotation. Substitut- 
ing the standard values, a = 8x 10~ 3 and C w = 2 x 10~ 5 , 
in the global model we have rd w = 1.4 AU. 

The disk winds in the inner region, r < r dw , do not 
have sufficient energy to escape from a disk by accre- 
tion. The fates of these wind materials are (i) accreting 
directly to a central star if the angular momentum is re- 
moved, (ii) blown away by the stellar winds (see §5.4.3), 
or (iii) returning back to the disk (see Figure 18 for the 
schematic picture). If the most of the wind gas follows 
the processes (i) or (ii), our calculations in §4 give the 
correct surface gas densities. On the other hand, if the 
process (iii) is dominant, we overestimate the escaping 
mass flux of the disk winds. We further discuss the es- 
cape of the disk winds in §5.4. 



time=0 , 




r(AU) 

Fig. 13. — Time evolution of disk surface density. The black thin 
lines are the results without disk winds (Model I); the blue thin 
lines are the results with the disk winds of the constant mass flux, 
C w = C Wj fl (Model II); the red thick lines the results with the 
disk winds which take into account the energetics limiter for C w 
(Equation 24). The dotted lines are the initial values. The dashed, 
solid, and dot-dashed lines are the results at 10 5 , 10 6 , and 10 7 yrs. 

Model II (constant C w ; blue thin lines). As expected, the 
difference between the red and blue lines appears only in 
r < rd w (= 1.4 AU), the surface densities in the outer 
region are the same. Although the surface density in the 
inner region becomes larger in the case with the C w lim- 
iter, it is still much lower than the surface density of the 
no wind case. Then, the disk winds still play an essential 
role in the dispersal of protoplanetary disks even though 
taking into account the energetics of the disk winds. 



5.2. Global Modeling with Energetics 



We can take into account the energetics argument 
(§5.1) in the global model. We adopt the following lim- 
iter for disk wind mass flux, C w : 



C w = min(C w ,n,C tt 



(24) 



where C Wj0ng corresponds to the mass flux that gives the 
left-hand side of Equation (21) equal to zero, namely 

C w ,( 

_18_ 

r dr 
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\dt 


2 



rfi^(r 2 £ac 2 ) + r 2 £ftac 2 
dr s 



r 2 n 2 p mid c 6 



(25) 

Equation (24) reduces the mass flux in r < rdw to the 
value available from the liberated gravitational energy by 
accretion. In r > ?'d w , Equation (24) gives C w — C w a. 

Figure 13 displays the results which take into account 
the limiter for C w (thick red lines) in Model II, in com- 
parison with Model I (no disk wind; black thin lines) and 



5.3. Mass Loss Rate 

Figure 14 shows the mass accretion rate, M r , and the 
mass loss rate of disk winds, M z , which are respectively 
defined as 



and 



M z (r) = 2tt 



— 27rrSw r 



r'dr'(pv z ) v 



(26) 



(27) 



In the case with disk winds, the mass accretion rate de- 
creases for decreasing r because the mass is lost by the 
disk winds. The total mass loss rate by the disk winds 
is M z {ri n ) = 2.5 x 10~ 9 A/ yr" 1 and the mass loss rate 
from r > r c \ w that can be accelerated to infinity by ac- 
cretion energy is M z (rdw) = 1-4 x 1CP 9 Ai© yr _1 ; more 
than the half of the total wind mass loss can escape from 
the disk by the liberated gravitational energy. The mass 
loss by the disk winds is larger than the accretion rate. 
Assuming the truncation of the disk at 5-10 stellar radii 
(~ 0.1 AU), M r w 2 x 10" 10 M Q yr" 1 , which can be re- 
garded as the mass accretion rate through magnetic flux 
tubes directly connecting to the central star (Kenyon, Yi, 
& Hartmann 1996). 
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Fig. 14. — Structure of mass loss/accretion rates of Model II in 
the main paper at t = 10 6 yr. The solid line is the mass loss by the 
disk winds (Equation 27) and the dashed line is the mass accretion 
rate (Equation 26) of the model taking into account the disk winds. 
The dotted line is the mass accretion rate of the model without disk 
winds for comparison. The arrow indicates the location of r<j w . 
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Fig. 15. — Time evolution of the mass loss rate of the disk wind 
(solid) and accretion rate (dashed) of Model II. Here, the wind 
mass loss rate is the integration of r > r^„, M z {r^), that can 
potentially escape from a central star gravity by accretion energy. 
The shown accretion rate is the maximum value at each time (e.g. 
for t = 10 6 yrs in Figure 14, the maximum accretion rate is ob- 
tained at r ft: 15 AU). 

The obtained mass loss rate by the MRI-driven disk 
winds is larger than the mass loss rate, < 10 -10 — 10~ 9 
Mq yr -1 , predicted by the UV photoevaporation (e.g. 
Matsuyama ct al.2003). Recently, photoevaporation by 
X-rays from a central star is also proposed. It is reported 



that the mass loss rate due to X-rays may be larger than 
that by the UV photoevaparation, where the calculated 
mass loss rates are rather uncertain (Ercolano, Clarke, 
& Drake 2009; Gorti, Dullemond, & Hollcnbach 2009; 
Owen et al. 2010). Significant difference of the disk wind 
mass loss from the photoevaporation processes by UV 
or X-rays is time evolution. The mass loss rate by the 
disk winds is correlated with the accretion rate (Figure 
15), because the energy source of the disk winds is the 
gravitational energy liberated by accretion. Therefore, 
the wind mass loss rate is larger at earlier times when 
the accretion rate is larger; the disk winds significantly 
contribute to the dispersal of the gas component of a 
disk from the beginning. As a result, an inner hole forms 
from the early epoch and its size is gradually expanding. 
On the other hand, the UV photoevaporation mechanism 
is significant at the later times after the sufficient mass 
dissipates by accretion. Therefore, an inner hole is antic- 
ipated at later time when the evaporating mass flux be- 
comes comparable to the mass accretion rate (Alexander 
et al. 2006). The X-ray photoevaporation, which may 
be effective from slightly earlier time, is also expected 
to give the similar trend to the UV photoevaporation 
(Gorti, Dullemond, & Hollenbach 2009). 

5.4. Escape of Disk Winds 

We further continue the discussions on the escape of 
disk winds. 

5.4.1. Local Simulations with Larger Vertical Boxes 

In order to study the acceleration of the disk winds at 
higher altitudes, we perform the local 3D MHD simula- 
tions with larger vertical boxes. We here use the realistic 
vertical gravity, 

GM+z 2 r 3 

® z ( r 2 _|_ z 2)3/2 Z (j.2 _|_ 2 2)3/2 ' ^ ' 

where r is radial location from a central star. In the 
usual local simulations, the only leading term, g z ~ £l 2 z, 
is considered, because r does not appear explicitly and 
it can be treated more easily (Hawley et al.1995). On 
the other hand, the escape velocity is not defined at the 
expense of the simplification, and the vertical gravity 

is overestimated by a factor of ^ 2+ ^ 2 ^ 3/2 , which makes 

density at a high altitude unrealistically lower. To avoid 
these shortcomings we use Equation (28) by explicitly 
setting r (Table 2). Note that r = 20H corresponds to 
the location of r ps 1 AU for the MMSN. In these runs, 
we set the net vertical fields with /3 z .,„id = 10 6 at the 
midplanes. 



Model 


r 


Box Size 


II (Reference) 




-AH < z < AH 


VI 


20 


-6H < z <6H 


VII 


20 


-8H < z <8H 


VIII 


10 


-12H < z < 12H 



TABLE 2 

Radial positions and vertical box sizes of the local 
simulations. 

Figure 16 compares the time averaged vertical struc- 
tures of these cases (Table 2). The top left panel shows 
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Fig. 16. — Comparison of the disk wind structures with the differ- 
ent box sizes. On the left, we display the time-averages of vertical 
velocity and density. On the right we compare the time-averages 
of magnetic energies and a values only in the — 6H < z < 6H 
region. The solid lines are the results with r = 20H and box 
sizes, -6H < z < 6H (Model VI) and -8H < z < 8H (Model 
VII). The dashed lines are the results with r = 10H and box size, 
-12H < z < 12H (Model VIII). The dotted lines are the results 
of the reference case (Model II; r — > oo and —AH < z < AH). On 
the right we use the thick lines for Model VII. 

that the onsets of the disk winds take place at higher 
altitudes in larger box cases; the results depend on the 
simulation box size. However, the mass flux of the disk 
winds (pv z ) do not show a monotonic behavior as pre- 
sented in the bottom panel of Figure 17. By increasing 
the box size from ±4H to ±6H, the mass flux decreases 
at first 7 . However, increasing the box size from ±6H 
to ±8H, the mass flux increases slightly. The largest 
box size case with the smaller gravity, r = 10H, shows 
larger mass flux than these smaller box cases. The mass 
flux of the disk winds is bound by a lower limit and dose 
not becomes further smaller even though we use a larger 
vertical box. 

The dependence of the mass flux on the simulation 
box size can be explained by the saturation of the am- 
plified magnetic fields (the top right panel of Figure 16). 
When we increase the box size from ±4H to ±6if, the 
escaping mass from the disk surfaces becomes smaller 
at first. Accordingly, the escaping magnetic flux of the 
toroidal (y) and radial (x) components decrease, because 

7 Although we do not take into account r in the case with the 
vertical box of —AH < z < AH, this effect is negligible for suffi- 
ciently small z. 
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Fig. 17. — The a values (top) and the mass flux of the disk winds 
(bottom) for the different sizes of the vertical simulation boxes. 
The triangles are the results of Model II (r — > oo), the squares are 
the results with r = 20 H (Models VI & VII), and the circles are 
the results with r = 10H (Model VIII). The dotted line is the level 
of C w , fl . 

the magnetic fields are frozen in the gas. The amplifica- 
tion of magnetic fields is balanced with the escape with 
the disk winds in addition to magnetic reconnections. In 
the larger box cases the escaping flux becomes smaller 
and the saturated level of magnetic fields increases. As a 
result of the larger magnetic pressure, the mass is lift up 
to higher locations (the bottom left panel of Figure 16). 
The increase of the density inhibits further decrease of 
the mass flux of the disk winds, pv z , in a self- regulated 
manner. 

In spite of the increase of the magnetic field strength 
for larger boxes, a z does not increase so much (the 
lower right panel of Figure 16). This is because the 
magnetic field of the large box cases around z m 1.5H 
(the locations of the peaks) is dominated by the co- 
herent toroidal component which does not contribute 
to the anisotropic Maxwell stress. Therefore, the inte- 
grated a(= J dzpa z / J dzp), which directly determines 
the global mass accretion rate, does not depend on the 
vertical box size (the top panel of Figure 17). 

The top left panel of Figure 16 shows that the average 
velocities of the disk winds do not still reach the escape 
speeds (= y/2rfl = V2 (jj) Hfl = 2 (§) c s ). However, 
our conservative choice of the floor values, C w ,a, probably 
gives a reasonable estimate (the dotted line of Figure 17), 
first because the mass flux seems bounded by the lower 
limit as explained above, and second because there are a 
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Fig. 18. — Schematic picture of dispersal of a protoplanetary disk 
by disk winds. The wind material in the outer region, r > rj w , can 
stream out solely by the gravitational energy by accretion. On the 
other hand, the wind material in the inner location cannot escape 
by itself from the gravity of a central star. A fraction of the winds 
directly accrete to the central star after lift up. A fraction of the 
disk winds may be accelerated by the dynamical pressure of the 
stellar winds. If the stellar wind flux is not strong enough, the disk 
wind material returns back to the disk after transported outward. 

couple of mechanisms that favor the escape of the disk 
winds but are not included in this paper (see below). 

5.4.2. Magnetocentrifugal Winds 

Our local simulations do not take into account the ac- 
celeration (momentum input) of the disk winds by cen- 
trifugal force with global magnetic fields. If poloidal 
(r — z components) magnetic field lines sufficiently in- 
cline with respect to an accretion disk, the gas can flow 
out along with the field lines by the centrifugal force 
(Blandford & Payne 1982; Kudoh & Shibata 1998). 
Such global magnetic fields will be common in protoplan- 
etary disks as a result of the contraction of cold molecular 
clouds with interstellar magnetic fields. Therefore, our 
local simulations probably underestimate the momentum 
input to the disk winds. 

5.4.3. Stellar Winds 

So far, we have focused on the disk winds driven by 
MHD turbulence. Central T Tauri stars are also ex- 
pected to drive stellar winds by the mass accretion from 
disks (Hirose et al. 1997; Matt & Pudritz 2005) and the 
surface convections (Cranmer 2008). Observations of T 
Tauri stars show that the outflow rates range from 10 -11 
to 10 -8 Mq yr -1 , which is typically 0.01 - 1 times of 
the mass accretion rates (White & Hillenbrand 2004). 
Although it is difficult to determine the exact launching 
locations of the observed outflows, some fractions are ex- 
pected to come from central stars (Edwards et al. 2006). 

At present the stellar wind is regarded to give a minor 
effect to the dispersal of the gas component of protoplan- 
etary disks (e.g. Shu, Johnstone, & Hollenback 1992) 
because the stellar winds almost slide along the disk 
surfaces (Matsuyama, Johnstone, & Hollenbach 2009). 
However, if disk materials are lifted up as we have shown 
so far, the ram pressure of stellar winds can directly push 
away the lifted up materials because of the large eleva- 
tion angle. 



We discuss the role of stellar winds from a simple mo- 
mentum balance. Let us consider a situation, in which 
the lifted up gas by disk winds floats above a disk and 
the radial stellar winds hit the floating gas. If we as- 
sume that both lifted up gas and stellar wind gas move 
together radially outward after the hitting, we can esti- 
mate the radial velocity of the moving gas, u cm b, from 
the momentum balance as 

(M z - m + WM+)v cmh = WM+v+, (29) 

where M z / m is the rate of the mass that is supplied 
from disk winds but float above a disk without suffi- 
cient energy, W is the fraction of the solid angle ob- 
scured by the lifted up gas by the disk winds, and M* 
and v*(~ 200 — 400 km s -1 ) are the mass loss rate and 
velocity of the stellar winds. As a typical example, we 
consider the result of Model II at t = 10 6 yr (Figure 

14). Then, M z<in = 1.1 x 10" 9 M Q yr -1 ; we assume 
the lifted up gas fill up to the height of r/2, which gives 

W = J^2~ 1(1/2) dcosd ~ °- 45 - Thc S as that is liftcd U P 
by the disk winds but does not have the sufficient energy 
will distribute in r < rd w (= 1-4 AU). Here, we compare 
Vcmb with the escape velocity at 1 AU, u e sc,o ~ 42 km 
s _1 as a typical condition. Substituting these values into 
Equation (29), if M* > 4 x lO" 10 Mq yr" 1 , v cmh > v csc , 
is satisfied and the lifted up material can be blown away 
by the stellar winds. 

If A/* is smaller than this value, the stellar winds can 
blow away the disk wind material at sufficiently high al- 
titudes where the density is low. The disk wind material 
at lower heights will move outward after hit by the stellar 
winds but again return back to thc disk without sufficient 
energy When the returning location becomes r > ?*d w , 
the gas finally flow out by the disk winds, (see Figure 18 
for the schematic picture). 



6. SUMMARY 

In this paper, we have shown that the MRI-driven pro- 
toplanetary disk winds disperse the gas component of 
disks from thc inside out. If net vertical magnetic fields 
with moderate strength exist, the disk winds and ac- 
cretion switch-on from the inner locations, which forms 
an expanding inner hole. This mechanism naturally 
explains observed transitional disks with inner holes. 
Model calculations that incorporate UV or X-ray pho- 
toevaporation with accretion also expects an inner hole 
at the later stage after the significant fraction of the gas 
disappears (Alexander, Clarke, & Pringle 2006; Gorti, 
Dullcmond, & Hollenbach 2009). The main difference 
of our mechanism from the photoevaporation processes 
is that the MRI-driven disk winds expect an inner hole 
from the early times and its size gradually grows from 
< 0.1 AU to several tens AU during the evolution of 
~ 10 7 years. 

Future high resolution observations by ALMA will be 
able to resolve inner holes with ~ a few AU at distance 
of 100 pc. We hope that observations of protoplanetary 
disks with various epochs will reveal the time-evolution 
of inner holes. 
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The dead zone does not affect the disk winds so much 
and the effect is oniy hmited at the eariy epoch (< 10 6 
years) of the disk evolution, because the MRI-drivcn disk 
winds arc driven from the surface regions with sufficient 
ionization degree. Even though a large dead zone forms 
around the midplane, the disk winds are also driven in- 
termittently with quasi-periodic cycles of 5-10 rotations 
as a result of the breakups of large-scale channel flows, 
similarly to the no dead zone simulation (SI09). The 
intermittency of the simulated disk winds should be ob- 
served as the variation of the disk surfaces, which might 
explain the observed large time variations of young stars 
(Wisnicwscki ct al. 2009; Muzerolle et al. 2009; Bary, 
Leisenring, & Skruskic 2009). 

The inside-out clearing of protoplanctary disks by the 
MRI-driven disk winds suppress the infall of boulders be- 
cause the outward force by gas pressure gradient is small. 
This is suitable condition for the formation of planetes- 
imals by aggregation, where we also need to examine 
sticking condition to study the actual growth of solid 



materials (e.g. Okuzumi 2009). The inside-out clearing 
may also increase the dust-to-gas ratio in the inner part 
of a disk, which is also favorable for the formation of 
planctcsimals by gravitational instability (Sckiya 1998; 
Youdin & Shu 2002; Johansen ct al. 2007). The mi- 
gration of planets is also suppressed in the inner region 
where the gas is dispersed at early times. Then, the in- 
ward migration of newly formed (proto-)planets stop at 
a certain location as shown in Figure 12. On the other 
hand, the gas remains in the outer region, so that the 
formation of gas planets can proceed there. 
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